package EP::EPC::Core::Data_Selection::exp;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::Error ':try';

use EP::EPC::Core::Clustering::exp;

use EP::Core::Data::Parser::exp::user;
use EP::Common::DBI::Dataset_Metadata_Type;

use PDL;
use PDL::NiceSlice;

use EP::Rclient;

# TODO: refactor subselection by row ids

# sub: process_selection
#   Entry point for selection processing.
# returns:
#   The new resulting dataset.
sub process_selection
    {
    my $self = shift;
    my ( $dataset, $query ) = @_;

    # Do something with the dataset in the mean time... This will modify the data associated with the dataset. On any
    # errors, should reset the data! (e.g., by delete'ing $dataset->{data}).
    process_subselection( $self, $dataset, $query );

    delete $dataset->{data} and throw EP::Common::Error ( -value => CUD_EMPTY_SELECTION )
        if $dataset->data->matrix->isempty;

    # Copy the dataset now. Note: this causes $sel_dataset to point to the same raw data as the original.
    my $sel_dataset = $dataset->copy( { name => undef } );
    $sel_dataset->data( $dataset->data() );

    EP::Core::Data::Parser::exp::user->_write_dataset_data( $sel_dataset, $sel_dataset->data );
    my $meta = $sel_dataset->data->process_metadata();

    # HACK. This is to force any future read of raw data on the original (pre-selection) dataset to re-read it from the
    # disk, rather than read the modified raw data from sel_dataset.
    delete $dataset->{data};

    $sel_dataset->data->raw_data->inplace->setnantobad();
    return ( $meta, $sel_dataset );
    }

# sub: process_subselection
#   Run the processing of different types of subselections
sub process_subselection
    {
    my $self = shift;
    my ( $dataset, $query ) = @_;

    subselect_gene_similarity
        ( $self
        , $dataset
        , $query->param('genes')
        , $query->param('N')
        , $query->param('distance')
        , $query->param('impute_row_avg')
        , $query->param('rows_regexp') . ""
        )
        if defined $query->param('genes') and $query->param('genes') ne '' and $query->param('N') > 0;

    subselect_column_headers
        ( $self
        , $dataset
        , $query->param('subselect_column_headers')
        , $query->param('cols_regexp') . ""
        , $query->param('cols_negate') . ""
        )
        if defined $query->param('subselect_column_headers') and $query->param('subselect_column_headers') ne '';

    subselect_row_headers
        ( $self
        , $dataset
        , $query->param('subselect_row_headers')
        , $query->param('rows_regexp') . ""
        , $query->param('rows_negate') . ""
        )
        if defined $query->param('subselect_row_headers') and $query->param('subselect_row_headers') ne '';

    subselect_na_values
        ( $self
        , $dataset
        , $query->param('subselect_na_values')
        )
        if defined $query->param('subselect_na_values') and $query->param('subselect_na_values') < 100;

    subselect_stdev
        ( $self
        , $dataset
        , $query->param('subselect_stdev')
        , $query->param('subselect_stdev_columns')
        , $query->param ('greatest_n_stdevs' )
        )
        if $query->param('subselect_stdev') > 0 or $query->param('greatest_n_stdevs') > 0;

    subselect_anova
        ( $self
        , $dataset
        , $query->param('factors') . ""
        , $query->param ('FDR')
        )
        if defined $query->param('subselect_anova');

    subselect_ebayes
        ( $self
        , $dataset
        , $query->param('factors') . ""
        , $query->param ('gene_number')
        )
        if defined $query->param('subselect_ebayes');
    }

# sub: subselect_gene_similarity
#   Subselects a given number of genes most similar to given one(s).
sub subselect_gene_similarity
    {
    my $self = shift;
    my ( $dataset, $genes, $N, $distance, $impute_row_avg, $regexp ) = @_;

    my $gene_ids = EP::EPC::Core::Clustering::exp::get_gene_ids( $dataset, $genes, $regexp );
    $dataset = EP::EPC::Core::Clustering::exp::create_similarity_search( $dataset, $distance, $gene_ids, $N, $impute_row_avg );

    my $matrix = $dataset->data->raw_data;

    open CLUST, "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $dataset->folder->location . "/$dataset->{clustering}"
        or die "Couldn't open clustering results file $dataset->{clustering}!";

    my %selected;
    while (<CLUST>)
        {
        next unless /CLUSTER (.*?) SIZE=(.*?)/;
        my ( $cluster_nr, $cluster_size ) = ( $1, $2 );

        # cluster members
        $_ = <CLUST>;
        my @these_selected_gene_ids = split /\s/;
        shift @these_selected_gene_ids;
        $_-- foreach @these_selected_gene_ids; # gene ids are 0-based

        @selected{@these_selected_gene_ids} = (1) x @these_selected_gene_ids;
        }
    close CLUST;

    my @selected_gene_ids = keys %selected;

    foreach my $row ( keys %{ $dataset->data->{row} } )
        {
        delete $dataset->data->{row}{$row} if 0 == grep { $row == $_ } @selected_gene_ids;
        $self->{session}->status ("Selecting rows (row $row)") if $row % 500 == 0;
        }

    my $selpdl = $matrix(:,long (sort { $a<=>$b} @selected_gene_ids) );

    warn $selpdl->info;
    $dataset->data->{matrix} = $selpdl;
    }

# sub: subselect_stdev
#   Removes rows when values do not fall outside of a certain stdev range in a certain percentage of the columns.
#   ... or just chooses genes with n greatest stdevs.
sub subselect_stdev
    {
    my $self = shift;
    my ( $dataset, $nr_sd, $pct_col, $n_stdevs ) = @_;

    my $matrix = $dataset->data->raw_data;

    my ( $means, $rms ) = (statsover $matrix)[0,1];
    my ( $selection, $c_selection );

    if ( $n_stdevs == 0 )
        {
        my $selpct = sumover( abs( $matrix - transpose($means) ) - transpose($rms) * $nr_sd > 0 ) / ( $matrix->dim(0) );
        ( $selection, $c_selection ) = which_both( $selpct >= $pct_col/100 );
        }
    else
        {
        # return if number of requested stdevs is more than number of elements 
        return if $n_stdevs > $rms->nelem;

        my $sorted_rms = qsorti($rms);
        $sorted_rms = $sorted_rms( -1:0 );
        ( $selection, $c_selection ) = ( $sorted_rms ( sequence($n_stdevs) ), $sorted_rms ( $n_stdevs + sequence( $rms->nelem - $n_stdevs ) ) );
        }

    delete $dataset->data->{row}{$_} foreach ( list $c_selection );

    my $selpdl = ( $selection->isempty ) ? null : $matrix( :, qsort $selection );
    $dataset->data->{matrix} = $selpdl;
    }

# sub: subselect_na_values
#   Remove rows above a certain percentage of missing values.
sub subselect_na_values
    {
    my $self = shift;
    my ( $dataset, $pct_na ) = @_;

    my $matrix = $dataset->data->raw_data;

    my $nbadrpct = nbadover($matrix) / $matrix->dim(0);
    my ( $selection, $c_selection ) = which_both( $nbadrpct < $pct_na/100 );

    delete $dataset->data->{row}{$_} foreach ( list $c_selection );

    my $selpdl = ( $selection->isempty ) ? null : $matrix( :, qsort $selection );
    $dataset->data->{matrix} = $selpdl;
    }

# sub: subselect_row_headers
#   Process regular expressions on row headers.
sub subselect_row_headers
    {
    my $self = shift;
    my ( $dataset, $match, $regexp, $negate ) = @_;
  
    warn "REGEXP: $regexp, NEGATE: $negate";

    my $matrix = $dataset->data->raw_data;

    $match = EP::EPC::Core::Data_Selection::define_regexp($match)
        if $regexp;
    $match = qr/$match/ if $regexp; # compile regexp only once - here

    my @matches = split /\r\n|\r|\n/, $match
        if !$regexp;
    warn "Matching against " . join ( "; ", @matches );

    my @selected_row_ids;
    foreach my $row ( keys %{ $dataset->data->{row} } )
        {
        my $full_annot = $dataset->data->{row}{$row}{full_annot};
        my @annot_vals = $full_annot->Values;

        my $is_row_matched = 0;

        if (!$regexp)
            {
            # compute intersection between requested matches and all annotation for row/col
            my %union; my %isect; my %union2; my %isect2;
            foreach my $e (@annot_vals) { $union{$e}++ && $isect{$e}++ }
            foreach my $e (keys %union, @matches) { $union2{$e}++ && $isect2{$e}++ }

            $is_row_matched = scalar keys %isect2;
            }
        else
            {
            # compute number of matches of the regexp to all annotation vals for row/col
            $is_row_matched = grep { $_ =~ /$match/ } @annot_vals;
            }

        $is_row_matched = (!$is_row_matched) if $negate;

        if ($is_row_matched)
            {
            push @selected_row_ids, $row
                if 0 == grep { $row == $_ } @selected_row_ids;
            }
        else
            {
            delete $dataset->data->{row}{$row};
            }

        $self->{session}->status( "Selecting rows (row $row)" ) if $row % 500 == 0;
        }

    my $selpdl = null;
    $selpdl = $matrix(:,long (sort { $a<=>$b} @selected_row_ids) ) if scalar @selected_row_ids > 0;

    warn "SELECTED " . $selpdl->info;
    $dataset->data->{matrix} = $selpdl;
    }

# sub: subselect_column_headers
#   Process regular expressions on column_headers.
sub subselect_column_headers
    {
    my $self = shift;
    my ( $dataset, $match, $regexp ) = @_;

    my $matrix = $dataset->data->raw_data;

    $match = EP::EPC::Core::Data_Selection::define_regexp($match) if $regexp;
    $match = qr/$match/ if $regexp; # compile regexp only once - here

    my @matches = split /\r\n|\r|\n/, $match if !$regexp;

    my @selected_col_ids;
    foreach my $column ( keys %{ $dataset->data->{column} } )
        {
        my $full_annot = $dataset->data->{column}{$column}{full_annot};
        my @annot_vals = $full_annot->Values;

        my $is_col_matched = 0;

        if (!$regexp)
            {
            # compute intersection between requested matches and all annotation for row/col
            my %union; my %isect;
            foreach my $e ( @annot_vals, @matches )
                {
                $union{$e}++ && $isect{$e}++;
                }

            $is_col_matched = scalar keys %isect;
            }
        else
            {
            # compute number of matches of the regexp to all annotation vals for row/col
            $is_col_matched = grep { $_ =~ /$match/ } @annot_vals;
            }

        if ($is_col_matched)
            {
            push @selected_col_ids, $column if 0 == grep { $column == $_ } @selected_col_ids;
            }
        else
            {
            delete $dataset->data->{column}{$column};
            }

        $self->{session}->status ("Selecting columns (column $column)");
        }

    my $selpdl = $matrix( long( sort { $a<=>$b } @selected_col_ids ), : );

    warn $selpdl->info if DEBUG;
    $dataset->data->{matrix} = $selpdl;
    }

# sub: subselect_anova
#   Filters for genes that show up as diff. exp. under an ANOVA analysis
# params:
#   dataset - EP::Common::DBI::Dataset to filter
#   factor - which factor to use
#   FDR    - false discovery rate
sub subselect_anova
    {
    my $self = shift;
    my ( $dataset, $factor, $FDR ) = @_;

    die "Please select a factor for running ANOVA!" if $factor eq "";

    my $dataset_binfile = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $dataset->folder->location . "/" . $dataset->filename . ".bin";
    $factor += 2; # in R the factors will be +2, because of 1-based indexing and the ID column considered as a factor, too.

    my @anova = EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/tmpl/R/ep.gene_filter.tmpl.R");
ep.anova( "$dataset_binfile", $factor, $FDR );
__RCODE__

    my $selidx = (pdl @anova) - 1;
    my $matrix = $dataset->data->raw_data;
    my $selpdl = $matrix(:,qsort $selidx);

    foreach my $row ( keys %{ $dataset->data->{row} } )
        {
        delete $dataset->data->{row}{$row} unless any $selidx == $row;
        }

    warn $selpdl->info;
    $dataset->data->{matrix} = $selpdl;
    }

sub subselect_ebayes
    {
    my $self = shift;
    my ( $dataset, $factor, $gene_number ) = @_;

    die "Please select a factor for running ebayes selection!" if $factor eq "";

    my $dataset_binfile = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $dataset->folder->location . "/" . $dataset->filename . ".bin";
    $factor += 2; # in R the factors will be +2, because of 1-based indexing and the ID column considered as a factor, too.

    my @ebayes = EP::Rclient->R(<<__RCODE__);
app_root = "$EP::Config::h->{ROOT}";
source("$EP::Config::h->{ROOT}/cgi-lib/tmpl/R/ep.gene_filter.tmpl.R");
ep.ebayes( "$dataset_binfile", $factor, $gene_number );
__RCODE__

    my $selidx = (pdl @ebayes) - 1;
    my $matrix = $dataset->data->raw_data;
    my $selpdl = $matrix(:, qsort $selidx);

    foreach my $row ( keys %{ $dataset->data->{row} } )
        {
        delete $dataset->data->{row}{$row} unless any $selidx == $row;
        }

    warn $selpdl->info;
    $dataset->data->{matrix} = $selpdl;
    }

1;
